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Abstract 

Background: To understand the ecological and evolutionary consequences of species interactions in food webs 
necessitates that interactions are properly identified. Genetic analyses suggest that many supposedly generalist 
parasitoid species should rather be defined as multiple species with a more narrow diet, reducing the probability 
that such species may mediate indirect interactions such as apparent competition among hosts. Recent studies 
showed that the parasitoid Asecodes lucens mediate apparent competition between two hosts, Galerucella tenella 
and G. calmahensis, affecting both interaction strengths and evolutionary feedbacks. The same parasitoid was also 
recorded from other species in the genus Galerucella, suggesting that similar indirect effects may also occur for 
other species pairs. 

Methods: To explore the possibility of such interactions, we sequenced mitochondrial and nuclear genetic markers 
to resolve the phylogeny of both host and parasitoid and to test the number of parasitoid species involved. We 
thus collected 139 Galerucella larvae from 8 host plant species and sequenced 31 adult beetle and 108 parasitoid 
individuals. 

Results: The analysis of the Galerucella data, that also included sequences from previous studies, verified the five 
species previously documented as reciprocally monophyletic, but the Bayesian species delimitation for A lucens 
suggested 3-4 cryptic taxa with a more specialised host use than previously suggested. The gene data analyzed 
under the multispecies coalescent model allowed us to reconstruct the species tree phylogeny for both host and 
parasitoid and we found a fully congruent coevolutionary pattern suggesting that parasitoid speciation followed 
upon host speciation. 

Conclusion: Using multilocus sequence data in a Bayesian species delimitation analysis we propose that 
hymenopteran parasitoids of the genus Asecodes that infest Galerucella larvae constitute at least three species with 
narrow diet breath. The evolution of parasitoid Asecodes and host Galerucella show a fully congruent coevolutionary 
pattern. This finding strengthens the hypothesis that the parasitoid in host search uses cues of the host rather than 
more general cues of both host and plant. 
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Background 

Consumer species feeding on multiple resources may in- 
directly mediate interactions among their resources [1-3], 
and thereby affect both trait evolution of interacting spe- 
cies and eventually speciation. A classic example is appar- 
ent competition, where one species affects the density of 
other species through an increased density of the shared 
consumer [4-9]. We are interested in the potential for 
such interactions among a set of chrysomelid beetles 
(Galerucella spp.) and their parasitoids, including the 
ecological and evolutionary consequences. Previous 
studies within this system show that parasitism rates 
may be very high, and that parasitoids may mediate ap- 
parent competition between hosts [10-12]. The litera- 
ture suggests that the same parasitoid species Asecodes 
lucens (Nees) (Hymenoptera: Eulophidae) attacks all 
Galerucella species [10,11,13,14], and parasitoid mediated 
interaction could therefore involve additional species. 
However, cryptic species are common in parasitoid sys- 
tems and diet breadth may be more narrow than literature 
information suggest [15-18]. For instance, Smith et al. [15] 
showed that 17 tachinid fly species that were thought to 
be broad generalist parasitoids on lepidopteran larvae con- 
sist of 32 species with a more narrow diet breadth. Simi- 
larly, Apanteles leucostigrnus, thought to be a generalist 
parasitoid on at least 32 species of Hesperiidae caterpillars, 
revealed 36 distinct Barcode clusters interpreted as 
provisional specialist species by Smith et al. [19]. 

DNA-based species delimitation provides a novel 
method for discovering cryptic species where supposedly 
generalist consumers are shown to be morphologically 
more or less identical specialist consumers. However, 
the standard DNA barcoding protocols used to infer pu- 
tative species are based on single loci and do not test 
species limits with models tracking speciation and popu- 
lation genetic processes e.g., [20]. An improved protocol 
such as the General Mixed Yule Coalescent model (here- 
after GMYC) treats species-to-population transitions in 
a gene tree [21]. However, GMYC also assumes recipro- 
cal monophyly of species and that the gene tree is 
known without error, ignoring that non-monophyletic 
species in gene trees are common [22,23] and that the 
time needed for complete lineage sorting of ancestral 
polymorphism is substantial ([24], see also [25,26]). 
When deriving the species tree, one approach is to con- 
catenate multiple loci into a superlocus by assuming that 
gene trees are consistent among each other and with the 
species tree. An alternative approach is to explicitly 
model the coalescent process of each gene tree separ- 
ately under the constraint of a common species tree 
[27-31]. This later approach is implemented in BPP, a 
Bayesian species delimitation method [32], that in a re- 
cent evaluation came out as the most accurate species 
delimitation method under varying conditions of 



divergence time, number of loci and migration [33] . BPP 
offers an exciting new tool to detect recent speciation 
events and cryptic species and is particularly suitable for 
host-parasitoid systems due to the prerequisite of an 
informed prior (a "guide tree"; Yang and Rannala 2010, 
see material and methods). 

Cryptic species are thought to be more common in 
groups where chemical senses are more highly developed 
than vision, such as most insects, because changes in 
chemical communication do not necessitate morpho- 
logical changes for the evolution of reproductive barriers 
between sibling species [34]. Cryptic speciation in parasitic 
hymenoptera is particularly interesting as this group is 
often reported to mediate indirect interactions between 
host species [5,35-37], but also because host finding in 
parasitoids involves both host volatiles and volatiles from 
the plant associated to the herbivore host [38]. Host find- 
ing is often described along a detectability-reliability axis, 
where plant volatiles associated to the host are more 
abundant than host volatiles, but also have a lower reliabil- 
ity. Parasitoids using a combination of host and plant vola- 
tiles may then face greater difficulties in finding hosts on 
different host plants compared to parasitoid species that 
only use host volatiles for host finding. As a consequence, 
generalist parasitoid species may have to use more general 
cues for host finding and this may, in turn, result in less 
efficient host search cf. [39]. Due to such trade-offs, spe- 
cialisation and speciation cannot be separated from the 
chemical ecology underlying host finding. 

The chrysomelid beetle genus Galerucella contains a 
mixture of monophagous and polyphagous herbivore spe- 
cies that occur commonly in wetlands across the northern 
hemisphere. Within this genus, in contrast to the generally 
accepted dogma for speciation in herbivore insects, host 
shifts have not necessarily involved closely related host 
plant species. Speciation in Galerucella has instead oc- 
curred between co-occurring wetland plants of such 
distant relations as Salicaceae, Primulaceae, Lythraceae, 
Rosaceae, Polygonaceae, Nymphaeaceae and Betulaceae. It 
is not likely that such distantly related plant species share 
much volatile profiles beyond green leaf volatiles [40], and 
even these are likely to occur in very different proportions 
(supported by unpublished chemical analyses in our re- 
search group). One may then ask how parasitoids face this 
chemical diversity and yet manage to find hosts and do 
host shifts. Most Galerucella species are only parasitized 
during the larval stage by Asecodes spp, but some species 
are also attacked by other hymenopteran taxa. The female 
parasitoids mainly attack the early larval stages, and lay 
multiple eggs in one larvae. In our study area, hosts and 
parasitoids typically have one generation per year but 
exceptions do occur. 

In this paper, we use multilocus sequence data to test 
whether A. lucens attacking Galerucella larvae represent 
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one generalist or multiple specialist species. For this 
purpose, we collected Galerucella larvae from most of their 
respective host plant species, reared them to adults to col- 
lect the parasitoids and sequenced them. For comparative 
purposes, we also sequenced three genes of the collected 
Galerucella species, to complement previous analyses [41]. 
Bayesian species delimitation suggested that A. lucens con- 
sists of 3-4 cryptic taxa with a more specialised host use 
than previously suggested. The gene data analyzed under 
the multispecies coalescent model allowed us to construct 
the species tree for both host and parasitoid and we found 
a fully congruent revolutionary pattern suggesting that 
parasitoid speciation followed upon host speciation. 

Methods 

Study species and sampling 

Specimens of Galerucella (Coleoptera) and Asecodes 
lucens (Nees) (Hymenoptera) were collected in Sweden 
and Finland during 2011 (Figure 1, for details see 
Additional file 1: Table SI-1). Larvae of Galerucella were 
collected on their host plant, and subsequent rearings 
provided either adult beetle or parasitoid specimens in 
the lab. One species, G. nymphaeae, was collected but 
was not found to be parasitized and is therefore ex- 
cluded from the data set. All specimens were stored in 
99% ethanol. The sampling in Sweden was mainly from 
Uppland and Vastmanland but parasitoids infesting G. 
tenella and G. calmariensis were also sampled from 
northern localities. In particular, a systematic sampling 
was performed in Skeppsvik archipelago in Umea, where 
parasitism rates are typically much higher than in southern 
localities (> 70% vs. < 10%). 

As outgroups to the Galerucella dataset we downloaded 
sequences of Trirhabda bacharidis, Pyrrhalta spp. and 
Ophraella communa (Galerucinae, Chrysomelidae) from 
GenBank. Galerucella sequences available on GenBank 
were also downloaded and included in our dataset (see 
Additional file 1: Table SI-1). In the Asecodes dataset we 
included part of the COl sequence of Horismenus 
missouriensis, H. petiolatus and Pediobius spp. (Entedo- 
ninae, Eulophidae) as well as Quadrastichus haitiensis 
(Tetrastichinae, Eulophidae) as outgroup species (down- 
loaded from GenBank). Two Pediobius specimens were 
also reared from G. sagittariae, an host association that 
was not previously recorded (Additional file 1: Table SI-1). 

DNA extraction, amplification and sequencing 
Extractions 

Specimens for DNA extraction were chosen to represent 
samples from different collection sites and, in the case of 
G. sagittariae, to represent different host plant prefer- 
ences. We used QIAamp® DNA Mini Kit (QiaGen) for 
the DNA extraction, where one to two legs per specimen 
were used for Galerucella while two to four legs were used 



for Asecodes. Legs were put in ATL lysis buffer and pro- 
teinase K and left in 56°C overnight. The purification pro- 
cedure was carried out according to the manufacturers 
protocol except that the Asecodes extracts were eluted in 
100 ul AE buffer. 

DNA amplification 

For both Asecodes and Galerucella we amplified the 
mitochondrial DNA of cytochrome oxidase 1 (COl) 
gene and two nuclear genes; the internal transcribed 
spacer (ITS2) from the nuclear ribosomal region and the 
D2 region of the 28S ribosomal subunit. COl was ampli- 
fied in two fragments without overlap using the primer 
pairs LCO-1490 & HCO-2198 [42] and Cl-J-2183 (Jerry) 
& TL2-N-3014 (Pat) [43]. For the amplification of the 
nuclear genes we used the primer pairs ITS2f & ITS2r 
[44] and 28S_D2_F & 28S_D2_R [44] respectively. For 
Asecodes we also amplified the nuclear fragment of 
phosphogluconate dehydrogenase (PGD) using the primer 
pair PGD_hym_3F & PGD_hym_intRb (Malm unpubl). 
PCR reactions were performed with Ready-ToGo™ PCR 
beads (Amersham Biosciences) in a total reaction vol- 
ume of 25 ul only adding 1 ul each of forward- and 
reverse primers (10 uM), 2 ul template and dH 2 0. The 
cycling profile started with a 5 min denaturation step at 
95°C, followed by 38 cycles of 30 s at 95°C, annealing for 
30 s at 49°C (ITS2), 50°C (COl), 60°C (28S), extension 
for 1 min at 72°C, and a final extension step of 8 min at 
72°C. For amplification of PGD, 40 cycles were run, 45 s 
was used in the denaturation step and annealing was 
performed for 40 s at 54°C in each cycle. The results 
were visualized on an agarose gel and stained with 
GelRed™ (Life Technologies). The PCR products were 
cleaned using Exonuclease I and FastAP (Fermentas) and 
the sequencing reactions were performed with BigDye™ 
Terminator ver. 3.1 Cycle Sequencing Kit (Applied 
Biosystems). The same primers as for the PCR reactions 
were used. The sequencing reactions were then cleaned 
with DyeEx 96 kit (QIAGEN) and run on an ABI Prism 
3100 Genetic Analyzer (Applied Biosystems). Genetic 
fragments were sequenced in both forward and reverse di- 
rections. To align the ITS sequences we used the program 
Mafft [45] and the iterative refinement method which 
incorporates global pairwise alignment information. All 
other sequences were aligned and edited by eye in BioEdit 
v7.0.9 [46]. 

Phylogenetic analyses 

Galerucella terminals were assigned to species based on 
which host plant they were sampled as well as morpho- 
logical identification. The Asecodes terminals were assigned 
to which host they parasitized; e.g., Asecodes I G. lineola 
refers to an Asecodes specimen hatched from G. lineola. 
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Figure 1 Sampling locations for Galerucella and Asecodes specimens in Sweden and Finland. Circles refer to sites were Golerucella were 
sampled and triangles refer to the collection sites of Asecodes. 



We used MrModeltest2 v.2.3 [47] in conjunction with 
PAUP [48] to select the best model for the Bayesian ana- 
lyses, although the substitution rate matrix was not se- 
lected a priori (see below). For COl the 1 st , 2 nd and 3 rd 
codon positions were analyzed separately. Substitution 
model details for the Galerucella and Asecodes datasets 
are specified in Additional file 1: Table SI-2. The Bayesian 
phylogenetic analyses were performed with MrBayes 3.1 
[49,50] at the online Bioportal server (University of Oslo, 



Norway). The Galerucella dataset consisted of 53 termi- 
nals representing eleven Galerucella species and two 
outgroup taxa for the COl + nuclear gene dataset and 19 
terminals representing five Galerucella species and two 
outgroup genera for the nuclear gene dataset. The 
Asecodes dataset (COl + nuclear genes) included 107 taxa 
and six outgroup taxa. Each gene was analyzed separately 
(except the 28s for Asecodes) as well as combined, to 
evaluate the relative information content in each gene 
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fragment. In the datasets including more than one par- 
tition (COl) or including several genes, the substitu- 
tion model was set to "mixed", which implements 
reversible-jump MCMC across the entire space of 203 
reversible 4x4 nucleotide substitution models [51]. 
All model parameters except branch lengths and topology 
were unlinked and relative rates between partions were 
allowed with a rate multiplier. We made two separate runs 
completing 10 million generations, each with four incre- 
mentally heated chains (T=0.2) and where sampling was 
done every 1000 th generation. The first 2 500 trees from 
each run were discarded as burn-in and the remaining 
samples pooled before calculating the majority-rule con- 
sensus. We checked that the separate analyses had con- 
verged using the average deviation of split frequencies 
diagnostic (< 0.02 in all runs), and the potential scale 
reduction factor (close to 1.00 for all parameters). 

We used the multispecies coalescent model [30] as 
implemented in *BEAST ('Starbeast') to infer the species 
tree from multiple gene trees for both Galerucella and 
Asecodes. In the multispecies coalescent model all model 
parameters are unlinked across loci, including the top- 
ology parameter, which allows the gene trees to differ in 
topology while being constrained by one and the same 
species tree. The Galerucella dataset consisted of 18 ter- 
minals from five species and for one mitochondrial 
(COl) and two nuclear loci (28S and ITS). Substitution 
models were set according to selected model for each 
loci by MrModeltest [47]. We applied a strict clock 
model on branch lengths and calibrated the COl partition 
to 0.0177 substitutions per site per million years after a re- 
cent COl clock rate estimate of another Polyphaga beetle 
family [52]. 

The Asecodes dataset consisted of 38 terminals from four 
species (as identified by the Bayesian species delimitation 
method, see below) and for one mitochondrial (COl) and 
three nuclear loci (PGD, ITS and 28S). We applied a strict 
clock model on branch lengths and calibrated the species 
tree by applying a normal distribution prior (mean = 4.4, 
stdev = 0.7) on the root node according to the 95% HPD 
interval of the estimated host root node age. The purpose 
of this calibration was to see if, under the assumption of 
contemporary speciation at the root node, branching events 
after the root node predate or postdate respective 
branching event in the host tree. While the latter would be 
in agreement with (but does not ascertain) cospeciation, 
the former would falsify any such hypothesis. Note that the 
substitution rate calibration in the Galerucella analysis is 
applied to the COl gene tree whereas the root node cali- 
bration in the Asecodes analysis is applied to the species 
tree. The latter was enabled by hand editing of the xml file 
following McCormack et al. [53]. 

Both analyses were run four times independently for 
100 million MCMC generations and sampled every 1000 



generation. Tracer [54] was used to examine the conver- 
gence across runs and the ESS values of sampled param- 
eters. After 10% of each run was discarded as a burn-in, 
the remaining samples from all runs were pooled and 
the maximum clade credibility species tree was calcu- 
lated using the mean node heights. 

Species delimitation analyses 

To test if the Asecodes parasitoid on Galerucella consists 
of one generalist species parasitizing multiple host species, 
or several specialists each attacking a single host, we used 
a Bayesian species delimitation method described by Yang 
and Rannala [32]. The method as applied in the software 
BPP uses reversible-jump MCMC to sample different 
species delimitation models and estimate the posterior 
probability of each model. BPP accommodates both 
uncertainty in gene tree estimations by the Bayesian 
framework, and incomplete lineage sorting of species 
via the coalescent process model (both in contrast to 
the GMYC method below). This statistical power in- 
creases both with the addition of loci and individuals 
per species [55]. BPP needs a guide-tree as input for the 
analysis and both the number of terminals in the guide- 
tree and its topology can influence the result [56]. BPP 
only estimates the posterior probability of all possible 
ways to collapse the nodes in the guide-tree into fewer 
species. It does not test any alternative groupings of 
individuals into species that cannot be derived from col- 
lapsing a node in the guide-tree, nor does it test to split 
the dataset into more species than the given terminals 
in the guide-tree. The dependency on a guide-tree 
effectively limits the search space and was a practical 
necessity introduced by Yang and Rannala [32]. As such, 
BPP is not suitable for blind biodiversity assessment, 
but explicit a priori species hypotheses need to exist to 
avoid artificial divisions like the smallest geographical 
sampling localities [56]. In this respect, parasitoids have 
the advantage of having natural and biologically relevant 
units, the host species, which in our analysis serve as 
delimiters of the maximal number of terminal units 
(five). We also tested guide-trees with four and three 
terminals, which only confirmed the relevant subset of 
results from the five terminal guide-tree analysis and 
will not be discussed further. The guide-tree topology 
was taken from the *BEAST analysis, all nodes of which 
were supported by a posterior probability of 1.0 and 
hence we did not need to test the effect on the species 
delimitation of any alternative plausible topology as 
guide-tree (compare with 56). 

To calculate the likelihood of different species de- 
limitations, a model with two types of parameters is 
needed (the number of parameters depend on the number 
of species in the model): G = 4Ne[t and r. G is a product 
of the effective population size and the mutation rate and 
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t is the root age measured in expected number of substi- 
tutions per site to the tips. Disregarding variations in the 
mutation rate, a relatively large value of G means a large 
effective population size and a relatively large value of r 
specifies an ancient divergence. Each terminal species in 
the species delimitation model and each ancestral inter- 
node (ancestral species) have a G parameter whereas each 
split between species has a r parameter. In a Bayesian 
framework both parameters need a prior which is speci- 
fied by a gamma distribution in BPP. The gamma distribu- 
tion is defined by two numbers, a and /?, and the mean of 
the distribution is alf$. The mean should not be orders of 
magnitude away from the posterior of the parameter as es- 
timated from the data, otherwise the prior can cause 
flawed results [55]. As a diffuse prior we used a = 1 in all 
analysis and tested various combinations of ft for G and r 
around the estimated posterior for the parameters. The es- 
timated posterior estimates were G ~ 0.05 and r ~ 0.02 
(see Additional file 1: Table SI-3). All analyses were there- 
fore run in four combinations of prior gamma distribu- 
tions: 1: G: G(l,10), r: G(l,10), 2: G: G(l,10), r: G(l,100), 
3: G: G(l,100), r: G(l,10) and 4: G: G(l,100), r: G(l,100), 
representing all combinations where the mean is between 
0.1 and 0.01 and which includes the posterior estimates 
for both parameters. We also tested the effect of the num- 
ber of loci on the power of the method. To do that we se- 
quentially added loci in the order 1 locus (COl), 2 loci 
(COl+PGD), 3 loci (COl+PGD+ITS), 4 loci (COUPGD+ 
ITS+28S). It could be argued that the adjacent ITS and 
28S gene segments should be analyzed as linked, but be- 
cause of the complex pattern of recombination and con- 
certed evolution in the tandemly arrayed rDNA loci [57], 
we prefer to analyze them as separate loci. If in doubt 
about the independence between the two loci, the result 
from the 3 loci analyses can be consulted. Each analysis of 
number of loci and prior combination was run 4 times 
independently to check the convergence of results and 
variation of estimated posterior probability for species 
delimitation models. Each analysis consisted of 50000 
MCMC generations sampled every 5 th generation and 
discarding 10% as burn-in. The two different algorithms 
for the reversible-jump MCMC described by Yang and 
Rannala [32] and named "0" and "1" gave similar results 
and we therefore only report the result from using the 
algorithm 0. 

The GMYC analysis was performed in R statistical pack- 
age with the help of ape, gee, MASS, paran and splits pack- 
ages. The GMYC method takes as input an ultrametric 
gene tree which we derived from calculating the maximum 
clade credibility tree with mean node heights from the 
sampled COl gene trees in a BEAST analysis using a strict 
molecular clock and a codon-partitioned model according 
to results from MrModeltest. The method calculates the 
likelihood of a mixed Yule-coalescence model applied with 



a single threshold across the tree at every node. At each 
threshold, deeper nodes are modeled as speciation events 
according to a Yule model whereas each group of younger 
nodes is modeled separately according to the coalescent 
process model. The maximum likelihood solution generally 
identifies the 'kink' or point of increased diversification rate 
in a lineage-through-time plot of trees with multiple indi- 
viduals per species for multiple species. This node serves 
as a species delimitation point under the assumption of 
species monophyly and not permitting any speciation event 
to be younger than the deepest coalescence. The maximum 
likelihood of the GMYC model was tested with a likelihood 
ratio test against a null model treating the entire tree as a 
single coalescent (i.e. against a one-species model). 

Results 

In total the Galerucella alignment consisted of 2766 
basepairs whereas the alignment of Asecodes contained 
3089 basepairs (see Additional file 1: Table SI-4 for 
details). 

Phylogenetic analysis 

The multispecies coalescent model for Galerucella 
gave a well supported species tree with the topology 
(G. sagittariae (G. lineola (G. tenella (G. pusilla + G. 
calmariensis)))). Subgenus Neogalerucella (including 
G. lineola, G. tenella, G. calmariensis and G. pusilla) 
and subgenus Galerucella (including G. sagittariae) 
were monophyletic as defined by Borghuis et al (2009), 
and G. lineola was basal in the Neogalerucella clade. All 
nodes had a posterior probability support of >0.99 and the 
age of the root node was estimated to 4.3 my (95% HPD: 
3.0 - 5.9) with the mitochondrial clock calibration 
(Figure 2). The standard concatenated model with uncon- 
strained branch lengths (non-clock model aka time-free 
model) of all genes combined resolved the species tree as 
the coalescent model to the same topology with high 
branch support (>0.99), but is here rooted using outgroups 
and includes a few additional species in the G. sagittariae 
group (Figure 3a). Analyzed alone, the mitochondrial COl 
gene tree is rooted differently compared to the nuclear 
(Figure 3b) and combined tree (COl and nuclear genes), 
which justifies the multispecies coalescent approach where 
gene trees are allowed to differ. It is interesting to note 
that the closely related G. calmariensis and G. pusilla 
were recovered as reciprocally monophyletic in the 
analyses of the nuclear genes but not by our COl data 
or by the mitochondrial data in Borghuis et al. [41]. The 
divergence date in the species tree is estimated to 77 
000 years ago (95% HPD: 19000-148000). 

The multispecies coalescent model for Asecodes with 
species defined according to the BPP analysis resolved 
the species tree to {Asecodes I G. sagittariae {Asecodes I 
G. lineola {Asecodes I G. tenella + Asecodes I G. pusilla I 
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Figure 2 Topologies of Galerucella and Asecodes species, analysed with multispecies coalescent models in *BEAST. The analyses are 
based on sequences of C01, 28S, ITS {Galerucella and Asecodes) and PGD {Asecodes). Branch labels indicate posterior probability values and the 
bars refer to the confidence interval of the node ages. The Galerucella dataset was rooted using a strict clock while the rooting of Asecodes was 
calibrated after the estimated age of Galerucella. The Asecodes specimens are coded as "Asecodes /the name of the Galerucella species that 
is parasitized". 



calmariensis)) with all posterior probability branch 
support 1.0 (Figure 2). Note that, as in the Galerucella 
analysis, the root was inferred from the clock model 
and not outgroups. As this is in perfect agreement with 
the phylogeny of the hosts, we calibrated the root node 
after the estimated age on the equivalent host chrono- 
gram node. For this we used the 95% highest posterior 
density of the host root node height to define a normal 
distribution (4.3, 0.7) as a prior on the parasite root 
node that takes the calibration uncertainty into ac- 
count. This resulted in the node for the split of 
Asecodes on G. lineola and the node for the split of 
Asecodes on G. tenella to postdate the equivalent nodes in 
the host tree. The opposite would have been a potential 
falsifier of a cospeciation scenario. 

With a standard concatenated model time-free (non- 
clock) model, and using outgroups to root the tree in a 
Bayesian analysis, three distinct clades were recovered cor- 
responding to parasitoids on G. sagittariae, G. lineola and 
G. tenella-calmariensis-pusilla. The Asecodes specimens 
that infest Galerucella sagittariae were recovered as sister 
to the remaining Asecodes with high posterior probability, 
in agreement with the multispecies coalescent model. 



Asecodes on G. sagittariae and on G. lineola were each re- 
covered as reciprocally monophyletic with high posterior 
probability support (1.0 and 0.91 respectively). Asecodes 
on G. tenella, G. calmariensis and G. pusilla however were 
not recovered as reciprocally monophyletic but were 
mixed in one large cluster, with strong branch support 
(1.0) (Figure 4). The COl gene tree was identical in top- 
ology to the combined tree whereas the nuclear dataset 
alone gave a completely unresolved topology in the 
ingroup because of few variable sites (Additional file 1: 
Figure SI-1). This gives an indication of the relative im- 
portance of different genes to the combined result. 28S 
had only a single variable site that separated Asecodes on 
G. sagittariae and G. lineola from the remaining parasit- 
oids. Two variable positions defined the G. sagittariae 
parasitoids in the ITS gene fragment. 

Species delimitation 

The BPP analyses yielded conclusive results and mul- 
tiple runs with the same settings were in close agree- 
ment. Independent of the prior combinations used, the 
species delimitation model with highest posterior 
probability was the four-species model (Figure 5). The 
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Figure 3 Phylograms from Bayesian analysis of Galerucella species, (a) Phylogram based on concatenated C01, 28S and ITS sequences 
(2677 basepairs). (b) Phylogram based on the nuclear sequences of 28S and ITS (1203 basepairs). Branch labels indicate posterior 
probability values. 



four-species model delimits a separate specialist 
Asecodes species each for G. sagittariae, G. lineola and 
G. tenella but a common species for G. pusilla and G. 
calmariensis. This model ranged in posterior probabil- 
ity from 0.73 to 0.99 (0.83 to 0.99 when analyzing all 
four genes) (Figure 5). The second best model was the 
five-species model where also Asecodes on G. pusilla 
and G. calmariensis was divided into separate species 
with a posterior probability of 0.002 to 0.27. The three- 
species model separating only Asecodes on G. sagittariae, 
G. lineola and remaining Galerucella had a marginal but 
non-zero support and the posterior probability for the 
two- and one-species model was 0 across all prior space 
investigated. 



The effect of the number of loci was most noticeable 
between one (COl) and two (COl + PGD) loci which 
yielded increased support for the four-species model. 
The third (ITS) and fourth (28 S) loci contained only 
three and one variable sites respectively and had there- 
fore marginal discriminatory power between species 
delimitation models. In fact, with the prior on 0 set to a 
relatively large population size G(l,10) the addition of 
ITS and 28S resulted in a decreased probability for the 
same model (Figure 5a). 

The prior on the root age x had very little effect on 
the posterior of delimitation models. With four loci and 
the prior combination 0: G(l,10), x: G(l,10) support 
for the best model was 0.83-0.84 and remained unchanged 
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Figure 4 Phylogram from Bayesian analysis of Asecode specimens, (a) Phylogram based on concatenated C01, PGD, 28S and ITS sequences 
(3124 basepairs). Due to short branches and lack of structure within the three clades, they are visualized as triangles codes as "Asecodes /the 
Galerucella species which are parasitized". The three clusters include 31 {"Asecodes/G. sagittariae"), 57 "Asecodes/G. tenella-pusilla-calmariensis" and 
19 "Asecodes/G. lineola" individuals, (b) Result of the GMYC analysis on an ultrametric C01 gene tree from BEAST. Splits of thick black branches 
indicate speciation events whereas splits of thin black branches indicate within-species coalescence events. Grey branches are ambiguous in the 
GMYC analysis and all solutions from 3 to 7 species are included in a +/- 0.5 Log likelihood confidence interval. The single coalescence 
one-species model is rejected in a likelihood ratio test p=0.011. 



with 0: G(l,10), t: G(l,100). The effect of the prior on 
0 was more noticeable and with the same priors on t 
but with G(l,100) on 0 instead of G(l,10), the posterior 
probability increased to 0.97-0.98 for the best model. 
The support for three separate Asecodes species, on 
G. sagittariae, on G. lineola and on the remaining 
Galerucella species, is conclusive and never received a 
posterior probability of less than 1.0. The support for a 
fourth species Asecodes on G. tenella was strong but not 
entirely conclusive, as the three-species model treating 
Asecodes on G. tenella, G. calmariensis and G. pusilla as 
a single species also received some support (up to a pos- 
terior probability of 0.1 in one analysis). In addition, as 
the interpretation of allopatric populations with limited 
gene flow is a controversial issue we tested the BPP 
method on a sympatric subset of the data with only all 
Asecodes on G. tenella and G. calmariensis sampled 
from the same localities in northern Sweden. Used on 
sympatric populations, species defined unambiguously in 
the BPP analysis should be uncontroversial [55]. For the 
same four prior combinations, the posterior probability 
was more evenly spread across the three-, four- and five- 
species models with zero probability for the one- and 
two-species models (Figure 5c). The five-species model 
received mean posterior support of 0.41-0.72. That the 
dataset consist of at least three species, as opposed to 
one, is therefore conclusive, but to be conservative 



against oversplitting, the three-species model cannot be 
rejected even though it is the least supported of the 
three. 

The GMYC model delimited five species with the 
threshold at the maximum likelihood solution (Figure 4b). 
The log likelihood of the GMYC model at the optimal 
threshold (247.2) was significantly better than the null 
model of a single coalescent (241.7) in a likelihood ratio 
test (p<0.011) and a one-species model could be rejected. 
However, the delimitation of three, four, five, six and seven 
species was less than 0.5 log likelihood units away from 
the maximum likelihood solution and should therefore be 
included within a reasonable confidence limit. In contrast 
to the BPP, the GMYC method does not accommodate in- 
complete lineage sorting, but assumes species monophyly. 
As the COl gene tree does not resolve Asecodes individ- 
uals parasitizing G. tenella, G. calmariensis or G. pusilla 
as reciprocally monophyletic, the delimitation of four or 
five species is not exactly the same as that assumed and 
tested in the BPP analysis. The GMYC method applied to 
the reduced sympatric dataset delimited three species at 
the maximum likelihood solution (data not shown). 

Discussion 

The analysis of mitochondrial and nuclear genetic markers 
indicate that Asecodes lucens is not one generalist species 
but at least three species with a more narrow host use. 
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Figure 5 Mean posterior probability value of species delimitation models. The priors @ and t refer to the ancestral population size and 
age of root respectively. The bars indicate standard deviation, a. A four-species model and different number of loci for priors on 0 and t. Loci: 
1 - C01; 2 - C01, PGD, 3 - C01, PGD, ITS; 4 - C01, PGD, ITS, 28S. b. Analyses using all four genes resulted in highest support for the four-species 
model followed by the five-species model, c. Analyses using all four genes but only including the Asecodes specimens infesting G. tenella and G. 
colmoriensis from the sympatric population in the northern part of Sweden (but including all other specimen infesting G. sogittorio, G. lineolo and 
G. pusilla). 



Two parasitoid species seem to be specialist on their re- 
spective host, Galerucella sagittariae and G. lineola. The 
three remaining Galerucella species are seemingly 
attacked by the same parasitoid species, even though these 
parasitoids also show some population differentiation. 
There were tendencies that the parasitoids attacking G. 
tenella, which feeds on F. ulmaria, had diverged from the 
parasitoids attacking G. calmariensis and G pusilla, which 
both feed on L. salicaria. The mitochondrial dating fur- 
thermore show that speciation in parasitoids are sequen- 
tial events following speciation in their host, and a recent 
but perhaps not completed split between parasitoids 
attacking G tenella and G calmariensis/pusilla, would 
follow this pattern. 

The derived phylogeny of Galerucella matches the 
previous analysis by Borghuis et al. [41] where the split 
between what is considered two subgenera, Galerucella 



and Neogalerucella, is most ancient. In our data set, the 
subgenus Galerucella only includes G. sagittariae 
whereas the other species belong to the subgenus 
Neogalerucella. Borghuis et al [41] based their analysis 
only on mitochondrial gene fragments (COl and 
NADH-2), and found strong support for the monophyly 
of G. tenella and G. lineola but could not resolve G. 
pusilla and G. calmariensis as reciprocally monophy- 
letic, indicating either a recent divergence or mitochon- 
drial introgression. Our study also included two nuclear 
DNA fragments (the D2 region in 28S and ITS2), which 
provides a test of the alternative explanations. Although 
our nuclear DNA datasets were small, divergence 
between G. pusilla and G. calmariensis was strongly 
supported. These two taxa both use Lythrum salicaria 
(Lythraceae) as host plant, and both ecological and 
morphological information suggest that the two species 



Hamback et al. BMC Evolutionary Biology 201 3, 13:92 
http://www.biomedcentral.com/1471-2148/13/92 



Page 11 of 14 



are reproductively isolated. The male copulatory organs 
are distinctly different, the body size differs and both 
larval and adult colour differs [58,59]. In these characters, 
there is quite small overlap. In addition, Nokkala and 
Nokkala [60] found karyotypic differences which further 
support species status. In contrast to the conclusion 
drawn by Borghuis et al. [41], based on the results from 
mitochondrial gene fragments, the lack of monophyly in 
the mitochondrial genes for G. pusilla and G. calmarensis 
might rather indicate some recent "phenomenon" such as 
introgression, Wolbachia infestation etc. However, a larger 
sample from each population is necessary to establish if 
lineage sorting is indeed complete for the nuclear markers, 
which theoretically should sort slower than mitochondrial 
counterparts. Due to the maternal inheritage of the mito- 
chondrion, effective population size is lower (l/4 th that of 
nuclear genes) and lineage sorting is therefore faster than 
for nuclear genes. 

The Bayesian species delimitation analysis for the A. 
lucens group was performed on the same set of nuclear 
and mitochondrial genes as the analysis for Galerucella, 
with the addition of the nuclear gene PGD, and provided 
strong evidence for population differentiation. It is pos- 
sible to conclude that A. lucens should be split at least 
into three species but additional data may strengthen 
the indication of further splits; the species delimitation 
analysis suggested 3-5 species. The molecular data have 
been confirmed by morphological studies that found dif- 
ferences in wing patterns among at least 3 taxa [61]. 
Comparisons with type specimens of available (synony- 
mized) names suggested the identity of two taxa, A. lucens 
parasitizing G. sagittariae and A. parviclava (Thomson) 
parasitizing G. tenella, G. calmariensis and G. pusilla, A 
third species, parasitizing G. lineola, represents a new 
species named A. lineophagum Hansson & Hamback [61]. 
There were no morphological characters supporting a 
further subdivision of A. parviclava. 

The BPP method distinguishes populations as different 
species if the per generation migration rate Nm « 1 [55], 
and the interpretation is therefore unambiguous in sym- 
patry but distance-decay patterns in allopatry may alone 
contribute to partial genetic isolation. Allopatry could ex- 
plain why the analysis with all data included identified 
population differentiation between parasitoids collected on 
G. tenella versus G. pusilla/ calmariensis. In the southern 
area, populations of G. tenella and G. pusilla/ calmariensis 
are often found in slightly different habitats and were 
typically collected in different localities. In the northern 
area, where population differentiation was seemingly 
weaker, G. tenella and G. calmariensis co-occur on Baltic 
shore lines and were collected in the same localities. It is 
possible that the southern population of A. parviclava is 
in an early stage of speciation, but the current data are 
insufficient to confirm this suspicion. 



The possibility of geographic variation in population dif- 
ferentiation of parasitoids is very interesting considering 
the previously documented parasitoid mediated indirect 
interactions between G. tenella and G. calmariensis in 
northern localities. Hamback et al. [11] found that para- 
sitism rates on G. tenella were higher at sites where this 
species co-occur with G. calmariensis compared with 
sites without G. calmariensis. The differences in parasit- 
ism rates translated into differences in the strength of 
interactions between the herbivore and its host plant 
[12]. In sites with L. salicaria and G. calmariensis , the 
attack rates on F. ulmaria by G. tenella were lower and 
the seed set were higher compared to sites without L. 
salicaria. Other studies show that this apparent compe- 
tition between the Galerucella species may have evolu- 
tionary consequences for F. ulmaria. Galerucella tenella 
has a very strong impact on plant fitness and studies on 
differently aged populations suggest that G. tenella 
causes a shift in the population towards plants with a 
lower height and with higher concentrations of potential 
antiherbivore compounds [62]. In the field, the lower 
quality of F. ulmaria as food for the larvae of G. tenella 
caused the beetle to expand its diet towards other 
Rosaceae plants [63] . Preliminary data suggest that these 
evolutionary changes in F. ulmaria only occur on 
islands without L. salicaria (Ericson & Stenberg, unpub- 
lished data). The hypothesis for these effects was that 
the parasitoids use both G. tenella and G. calmariensis, 
even though the behavioural mechanisms are not fully 
understood [64]. It was therefore important to know 
whether parasitoids collected from the two host were 
indeed the same population and this seems to be the 
case. The potential for a larger population differenti- 
ation in southern localities suggest that a similar appar- 
ent competition is not likely in these localities, and the 
parasitism rates in these localities are also typically 
much lower (<10% vs >70%). 

Speciation in the parasitoids for this system seem to 
follow the identity of the beetle host rather than the 
identity of the beetle s host plant. This finding suggests 
that host finding or recognition cues originate from the 
larval host rather than from the host plant, such as beetle 
produced pheromones or beetle specific plant cues. Two 
contrasts in particular supports this view. First, G. 
sagittariae feed on multiple host plant species that are not 
closely related, but we nevertheless found no population 
differentiation among parasitoids hatching from larvae 
collected on different plants. Second, G. sagittariae and G. 
tenella are known to feed on the same host plant, but we 
find that their respective parasitoid belong to different 
taxa. It seems less likely that plant produced volatiles pro- 
vide sufficient information for parasitoids to both differen- 
tiate between G. sagittariae and G. tenella, and at the 
same time to locate G. sagittariae on its different host 
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plants. Beetle produced compounds therefore seem more 
likely. Previous studies show that at least two Galerucella 
species (G. pusilla and G. calmariensis) produce aggrega- 
tion pheromones [65], but these compounds are produced 
by adult males and seem less likely to provide any in- 
formation on the whereabouts of Galerucella larvae. To 
further understand the speciation process in Asecodes 
parasitoids, we are currently working to identify the 
compounds involved during the search process. 

Speciation in parasitoids may follow different pathways 
according to published data, and the causes underlying 
this variability is not well understood. There are cases 
with host-association differentiation and seemingly tight 
cospeciation [66,67], similar to the one suggested in this 
study for the Galerucella- Asecodes system. In other 
cases, speciation in parasitoids is less well connected to 
the phylogeny of their host [68,69]. The current data 
both in this and most of the published cospeciation ex- 
amples cannot differentiate whether this is ecological 
speciation where hosts or parasitoid species are not geo- 
graphically isolated or cladogenesis where host and para- 
sitoids are isolated in pairs [70]. Differences in the 
speciation pattern may arise because the parasitism 
process involves several steps that create barriers to par- 
asitoids when their host switches diet. Diet switching 
and speciation in herbivorous insects is often connected 
to enemy free space where the host is more difficult to 
either find or exploit on an alternate plant [71-73]. Dif- 
ferent plants may produce quite different volatile profiles 
upon damage and this may cause problems for parasit- 
oid females to either find their host or even to identify 
the larvae as such. There are cases where a parasitoid 
species has geographic variability in the type of cues 
used during host search [74], but this area of research is 
poorly exploited. Besides affecting host search cues, 
plant quality is also known to affect herbivore immuno- 
competence and switching to an alternate plant may 
affect resistance to parasitoid attack. The Galerucella- 
Asecodes system seems ideal to study these processes in 
progress, both because host race formation is common 
and well documented among Galerucella beetles [75-77] 
and because interaction strengths between host and para- 
sitoid show quite large geographic differences. Further 
studies on species interactions in this system however ne- 
cessitates a better understanding also on the small scale 
population differentiation, to identify strengths of direct 
and indirect species interactions and their ecological and 
evolutionary consequences. 

Conclusion 

Our analysis, using multilocus sequence data in a Bayesian 
species delimitation analysis, confirms the phylogenetic 
structure within the chrysomelid beetle genus Galerucella 
but also that the parasitoids attacking Galerucella larvae 



should be subdivided into at least three species. This sub- 
division has later been confirmed with morphological data. 
The previous literature suggested that one species of the 
genus Asecodes attacked Galerucella larvae but our ana- 
lysis show that the parasitoids have a more narrow diet 
breadth. The evolution of parasitoid within Asecodes and 
the host Galerucella show a fully congruent revolution- 
ary pattern, even though the specific mechanism cannot 
be identified with the current data. The phylogeny of 
Asecodes show that host use is not connected to the plant 
species that host the beetle larvae as different host larvae 
on different plant species are attacked by different 
Asecodes species. This finding strengthens the hypothesis 
that the parasitoid uses cues of the host during search ra- 
ther than more general cues of both host and plant. Previ- 
ous studies within the Asecodes-Galerucella system has 
suggested that parasitoid-mediated interactions may be 
important for both beetle densities and plant performance. 
The observation that parasitoids in fact have a more nar- 
row diet breadth suggests that such indirect interactions 
may be restricted to subsets of the food web involving 
these two genera. 
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